(* Copyright 2007 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 *)
#light

namespace dnAnalytics.LinearAlgebra

open dnAnalytics.LinearAlgebra

/// A module which implements functional matrix operations.
module Matrix =
    
    /// Fold a function over all matrix elements.
    let fold (f: 'a -> float -> 'a) (acc0: 'a) (A: #Matrix) =
        let n = A.Rows
        let m = A.Columns
        let mutable acc = acc0
        for i=0 to n-1 do
            for j=0 to m-1 do
                acc <- f acc (A.Item(i,j))
        acc

    /// Fold a matrix by applying a given function to all matrix elements.
    let foldi (f: int -> int -> 'a -> float -> 'a) (acc0: 'a) (A: #Matrix) =
        let n = A.Rows
        let m = A.Columns
        let mutable acc = acc0
        for i=0 to n-1 do
            for j=0 to m-1 do
                acc <- f i j acc (A.Item(i,j))
        acc

    /// Create a 2D array from a matrix.
    let to_array2 (A: #Matrix) =
        let n = A.Rows
        let m = A.Columns
        Array2.init n m (fun i j -> (A.Item(i,j)))
    
    /// Checks whether a predicate holds for all elements of a matrix.  
    let forall (p: float -> bool) (A: #Matrix) =
        let mutable b = true
        let mutable i = 0
        let mutable j = 0
        while b && i < A.Rows do
            b <- b && (p (A.Item(i,j)))
            j <- j+1
            if j = A.Columns then i <- i+1; j <- 0
        b
    
    /// Chechks whether a predicate holds for at least one element of a matrix.
    let exists (p: float -> bool) (A: #Matrix) =
        let mutable b = false
        let mutable i = 0
        let mutable j = 0
        while not(b) && i < A.Rows do
            b <- b || (p (A.Item(i,j)))
            j <- j+1
            if j = A.Columns then i <- i+1; j <- 0
        b
    
    /// Checks whether a position dependent predicate holds for all elements of a matrix.
    let foralli (p: int -> int -> float -> bool) (A: #Matrix) =
        let mutable b = true
        let mutable i = 0
        let mutable j = 0
        while b && i < A.Rows do
            b <- b && (p i j (A.Item(i,j)))
            j <- j+1
            if j = A.Columns then i <- i+1; j <- 0
        b
    
    /// Checks whether a position dependent predicate holds for at least one element of a matrix.
    let existsi (p: int -> int -> float -> bool) (A: #Matrix) =
        let mutable b = false
        let mutable i = 0
        let mutable j = 0
        while not(b) && i < A.Rows do
            b <- b || (p i j (A.Item(i,j)))
            j <- j+1
            if j = A.Columns then i <- i+1; j <- 0
        b
    
    /// Map every matrix element using the given function.
    let map (f: float -> float) (A: #Matrix) =
        let N = A.Rows
        let M = A.Columns
        let C = A.Clone()
        for i=0 to N-1 do
            for j=0 to M-1 do
                C.[i,j] <- f (C.Item(i,j))
        C
    
    /// Map every matrix element using the given position dependent function.
    let mapi (f: int -> int -> float -> float) (A: #Matrix) =
        let N = A.Rows
        let M = A.Columns
        let C = A.Clone()
        for i=0 to N-1 do
            for j=0 to M-1 do
                C.[i,j] <- f i j (C.Item(i,j))
        C
    
    /// In-place assignment.
    let inplace_assign (f: int -> int -> float) (A: #Matrix) =
        for i=0 to A.Rows-1 do
            for j=0 to A.Columns-1 do
                A.Item(i,j) <- f i j
    
    /// In-place map of every matrix element using a position dependent function.
    let inplace_mapi (f: int -> int -> float -> float) (A: #Matrix) =
        for i=0 to A.Rows-1 do
            for j=0 to A.Columns-1 do
                A.Item(i,j) <- f i j (A.Item(i,j))
    
    /// In-place addition of two matrices.
    let inplace_add (A: #Matrix) (B: #Matrix) =
        for i=0 to A.Rows-1 do
            for j=0 to A.Columns-1 do
                A.Item(i,j) <- (A.Item(i,j)) + (B.Item(i,j))
    
    /// In-place subtraction of two matrices.
    let inplace_sub (A: #Matrix) (B: #Matrix) =
        for i=0 to A.Rows-1 do
            for j=0 to A.Columns-1 do
                A.Item(i,j) <- (A.Item(i,j)) - (B.Item(i,j))
    
    /// In-place componentwise multiplication of two matrices.
    let inplace_cptMul (A: #Matrix) (B: #Matrix) =
        for i=0 to A.Rows-1 do
            for j=0 to A.Columns-1 do
                A.Item(i,j) <- (A.Item(i,j)) * (B.Item(i,j))
    
    /// In-place scaling of every matrix element by a scalar.
    let inplace_scale (f: float) (A: #Matrix) =
        for i=0 to A.Rows-1 do
            for j=0 to A.Columns-1 do
                A.Item(i,j) <- f * (A.Item(i,j))
        
    /// Creates a sequence that iterates the non-zero entries in the matrix.
    let nonzero_entries (A: #Matrix) =
        seq { for i in 0 .. A.Rows-1 do
                for j in 0 .. A.Columns-1 do
                  if A.Item(i,j) <> 0.0 then yield (i,j, A.Item(i,j)) }
    
    /// Returns the sum of all elements of a matrix.
    let sum (A: #Matrix) =
        let mutable f = 0.0
        for i=0 to A.Rows-1 do
            for j=0 to A.Columns-1 do
                f <- f + A.Item(i,j)
        f
    
    /// Iterates over all elements of a matrix.
    let iter (f: float -> unit) (A: #Matrix) =
        for i=0 to A.Rows-1 do
            for j=0 to A.Columns-1 do
                f (A.Item(i,j))
        ()
    
    /// Iterates over all elements of a matrix using the element indices.
    let iteri (f: int -> int -> float -> unit) (A: #Matrix) =
        for i=0 to A.Rows-1 do
            for j=0 to A.Columns-1 do
                f i j (A.Item(i,j))
        ()
        
    /// Fold one column.
    let foldCol (f: 'a -> float -> 'a) acc (A: #Matrix) k =
        let mutable macc = acc
        for i=0 to A.Rows-1 do
            macc <- f macc (A.Item(i,k))
        macc
    
    /// Fold one row.
    let foldRow (f: 'a -> float -> 'a) acc (A: #Matrix) k =
        let mutable macc = acc
        for i=0 to A.Columns-1 do
            macc <- f macc (A.Item(k,i))
        macc

    /// Fold all columns into one row vector.
    let foldByCol (f: float -> float -> float) acc (A: #Matrix) =
        let v = new DenseVector(A.Columns)
        for k=0 to A.Columns-1 do
            let mutable macc = acc
            for i=0 to A.Rows-1 do
                macc <- f macc (A.Item(i,k))
            v.[k] <- macc
        v :> Vector
    
    /// Fold all rows into one column vector.
    let foldByRow (f: float -> float -> float) acc (A: #Matrix) =
        let v = new DenseVector(A.Rows)
        for k=0 to A.Rows-1 do
            let mutable macc = acc
            for i=0 to A.Columns-1 do
                macc <- f macc (A.Item(k,i))
            v.[k] <- macc
        v :> Vector